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Abstract 

The mechanisms underlying collective migration, or the coordinated movement of a population of cells, are not 
well understood despite its ubiquitous nature. As a means to investigate collective migration, we consider a wound 
healing scenario in which a population of cells fills in the empty space left from a scratch wound. Here we present 
a simplihed mathematical model that uses reaction-diffusion equations to model collective migration during wound 
healing with an emphasis on cell movement and its response to both cell signaling and cell-cell adhesion. We use 
the model to investigate the effect of the MAPK signaling cascade on cell-cell adhesion during wound healing after 
EGF treatment. Our results suggest that activation of the MAPK signaling cascade stimulates collective migration 
through increases in the pulling strength of leader cells. We further use the model to suggest that treating a cell 
population with EGF converts the time to wound closure (as function of wound area) from parabolic to linear. 
Keywords: Collective Migration, Cell signaling, Reaction-Diffusion Equations, Wound healing. 
Epithelial-to-Mesench5raial Transition 


1. Introduction 


Collective migration, or the coordinated movement of a population of cells is a ubiquitous biological process 
involved in tissue formation and repair, developmental regulation, and tumorigenesis Although there has been a 


heavy focus on collective migration in recent scientific literature, the mechanisms underlying this important process 


are still not well understood 




34| . Determining the underlying mechanisms behind collective migration 


behavior will likely provide insight into wound healing, as well as the epithelial-to-mesenchymal transition (EMT) 


B- 


and metastasis characteristics of tumor progression during cancer 

Significant intercellular communication is required for a group of cells to coordinate their movement as one unit 


protein kinase (MAPK) signaling cascade may promote collective migration 


affects many cellular activities, ranging from cellular migration to angiogenesis 



H 

1. The MAPK 
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IZ, 


481 . While several 


types of MAPKs appear to influence cell migration, the extracellular signal-regulated kinas 1/2 (ERK) MAPKs 
are dominant regulators of cellular migration. ERK influences cellular migration in a variety of cell types and is 
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downstream from the epidermal growth factor receptor (EGFR) |23l. 1271 l46l|. EGFR is a tyrosine kinase receptor 
that forms an active dimer state in the presence of ligand, and its main ligand is epidermal growth factor (EGF). 
Receptor dimerization initiates autophosphorylation of the intracellular EGFR tyrosine kinase domains, which 
recruit signal transducers, such as the Ras protein. These transducers then trigger various intracellular signaling 
cascades, most notably the MARK and ERK signaling cascades [l^ 3 |. Accordingly, we focus on the activation of 
ERK in response to EGF treatment in this study. 

Collective migration is defined by three main characteristics: physical and functional connection between cells. 


Q- 


multicellular polarity, and the deposition and remodeling of the extracellular matrix \m . We focus on the first 
of these hallmarks through cell-cell adhesion. We also note, however, that many aspects of single cell migration, 
such as molecular control of protrusions, polarity, shape generation, and cell-surface adhesion also affect collective 


migration [10l.ll3l.l2C 


M 


361 . Cell-cell adhesion is believed to have a stronger influence on collective migration than 


these aforementioned single cell features, yet its role is still poorly understood. [l2LIl3l|. For this reason, we ignore 
the effects of cell-surface adhesion and focus on cell-cell adhesion, as has been done in previous studies on collective 
migration fli Q,Q,Qi- 

Adherens junctions mediate cell-cell adhesion in tissues by attaching actin and intermediate filament cytoskele- 
tons between adjacent cells, thereby allowing for dynamic coordinate force transmission between the actin cytoskele- 
tons of physically contacting cells [l^ . Adherens junctions form between cells through transmembrane cadherin 
proteins, which bind actin filaments through zyxin, vinculin, and a-catenin adapter proteins [^, 261. Cadherin 
proteins have been reported as both promoters and inhibitors of cell migration, so their overall role is not clear 

innn 

12L 1 14 1181135| . Based on this uncertainty in the role of cell-cell junctions during collective migration and the 
MARK cascade’s well-established (stimulating) role js, 21, 31, 34|, we pose the following question: how does MARK 
activity influence the main role of cell-cell junctions during collective migration? 

To provide insight into collective migration, we will concentrate on wound healing induced in a sheet of cells 
from the human keratinocyte (HaCaT) cell line. Experiments, as previously described in j^, are conducted in 
which a wound comprised of open space is formed in a monolayer of cells by scratching away a portion of the 
population with a pipette tip. The remaining cell population begins a collective migration process into the wound 
and eventually fills in and closes the wound. Studies of this sort tend to split cells in the population into two 
distinct classifications: leader cells (those located at the wound interface) and follower cells (those located farther 
away from the wound). We use live-cell imaging to capture the cell profiles over time (every half hour) for 30 hours 
and analyze both untreated (denoted as mock) and EGF-treated cell populations to analyze the effects of both EGF 
and the MARK signaling cascade on collective migration. 

Mathematical models can provide insight into how EGF treatment stimulates collective migration, which may 
be difficult to decipher through experimental procedures. Some authors strongly advocate for mathematical models 
as necessary for a comprehensive understanding of signaling cascades due to their complex quantitative nature 
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. Accordingly, we focus on developing a simple first mathematical model to investigate the relationship between 
collective migration, the MAPK cascade, and cell-cell adhesion. In particular, we will develop two models to 
investigate how the MAPK signaling cascade affects cell-cell adhesion in response to EGF treatment during wound 
healing. 

Other studies have used various quantitative approaches to model collective migration. A partial differential 
equation (PDE)-based model is advantageous for modeling collective migration, as it allows for an in-depth analysis 
of the biol(^ical system paired with the ability to search for certain patterns, such as traveling wave solutions to 
the model ja, Il3|- This class of wave solutions frequently arise in models of biological systems 


3C 


3^,129 
3- Fis 


39| as well as 


Fisher’s equation 


in collective migration models due their ability to capture invasive wave front dynamicqj 

(originally used to model the spread of advantageous genes [3) is the classical example of a PDE model with 
traveling wave solutions, and it has been used previously to model collective migration in response to the wounding 
of a human peritoneal mesothelial cell population [30 1. 


A more recent study developed a two-dimensional continuum model based on the force of lamellipodia, cell- 
substrate adhesion, and cell-cell adhesion that can reproduce results from both contraction experiments of ente- 
rocytes and expansion experiments of Madin-Darby canine kidney (MDCK) cells . A later study demonstrated 
that simulations from this model also matched the leading edge velocity of experimental data in intestinal entero- 
cyte (lEC) cells [ 2 I. This model was further used as a means to identify the effects of wound shape and area on 
the time for wound closure. Ultimately, the authors concluded that the rate of change of the wound area during 
wound healing is likely proportional to the square root or first power of the wound area in lEC cells. This was an 
improvement over previous studies that assumed the rate of change of the wound area is constant, which yielded 
inaccurate estimates for larger wounds [3 - 


In 


38l |. the authors computationally recreated the results of 


34 


, i n which two waves of MAPK activation 
381 developed a four-species PDE model for 


were observed in MDCK cells during wound healing. The authors of 
MAPK activation during wound healing that can quantitatively recapitulate these two MAPK activation waves. 
The components of this model include an arbitrary diffusible ligand (such as EGF), its cell-surface receptor, an 
intracellular precursor to this ligand, and reactive oxygen species. This study ignored cell migration, however, and 
thus did not investigate the effect of these waves on wound healing. Other previous PDE models have focused on 
how cell signaling effects cell migration. In [3) the authors showed how one simple regulatory chemical (impacting 
only cell proliferation) is sufficient for a simple reaction-diffusion model to match data on wound healing of epidermal 
cells. 

In this study, we will develop two testable PDE models using reaction-diffusion equations in order to elucidate 
properties of the interaction between MAPK signaling, cellular migration, and cell-cell adhesion in a wounded 


^Discrete random walk models have also been proposed to model collective migration and can yield similar mean-field results to PDE 
models niia. 
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HaCaT cell population in response to EGF treatment. These models will aid in making predictions on cellular 
migration in epithelial cell populations. 

We use these two models to test two competing hypotheses regarding the influence of MAPK activation on the 
role of cell-cell adhesion during wound healing. The first hypothesis states that MAPK activity stimulates collective 
migration through decreases in the drag strength of follower cells. We denote the corresponding mathematical 
model as Model H since cell-cell adhesion will hinder migration. The second hypothesis states that MAPK activity 
stimulates collective migration through increases in the pulling strength of leader cells. We denote the corresponding 
mathematical model as Model P since cell-cell adhesion will promote migration. While both models are able to 
fit leading edge data, we demonstrate that the model based on the second hypothesis can better match various 
characteristics of the cell population during wound healing. These model results suggest that the dominant role 
of cell-cell adhesion in response to MAPK activation after EGF treatment is to promote leader cell pulling rather 
than inhibit the drag of follower cells. 

In Section[51 we present our two model derivations based on the above hypotheses. In Section[3J we demonstrate 
how Model P can match the leading edge dynamics better than Model H and then use Model P to predict wound 
closure as a function of wound area in HaGaT cell populations. Lastly, in Section [4l we conclude with a discussion 
on the implications of these results as well as plans for future work. 

2. Model Development 

Our mathematical models presented here consist of two coupled variables describing cell density, u(t, x), and the 
average cellular MAPK activation level, m(t), where our independent variables include time (t) and spatial location 
(cc). While we are investigating cell-cell interactions in this study, current techniques cannot directly measure these 
interactions. Hence, we use cell density as a means to investigate cell-cell interactions during migration, because it 
is measurable and indicative of physical interaction, so high cell densities should display a larger degree of cell-cell 
interaction effects on migration. We include the diffusion of cells and focus on its response to different levels of 
MAPK activity. This can be quantitatively described using a conservation law framework as: 

ut{t,x) =V ■ if^Vu), mt{t) = f^, (2.1) 

where the subscript t denotes differentiation with respect to time, fff denotes the rate of cellular diffusion and 

describes the activation rate of m[t). It should be noted that cell proliferation and death could be included in 
Equation (EH), but are assumed negligible over the course of the experimeniy. Decreases in MAPK activation could 
also be considered in Equation 1)2.II) . however, we do not see any noticeable decreases in MAPK levels during EGF 
experiments, so we exclude them as well. Note that while we simulate our model in only one spatial dimension due 

^The inclusion of a growth term had a minimal effect on model simulations and hypothesis evaluation (results not presented in this 
work). 
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to the one-dimensional plane wave-behavior in the wound healing data, the model can be easily extended to two or 
three dimensions. 

We will develop two different models in Section 12.11 Model H focuses on how cell movement towards the wound 
may be hindered by the drag of follower cells, while Model P focuses on how leader cells may pull on cells behind 
them to promote movement towards the wound. While we generally think of follower and leader cells in terms of 
their location in the population, we assume that each cell behind the first row in our population will either pull 
on those behind it or hinder those in front of it. The resulting models thus distinguish the dominant interaction 
(dragging or pulling) betweens neighboring cells in the population during wound healing. 

In the rest of this section, we describe the development of two simple models for different hypotheses on the 
emergence of collective migration during wound healing. In Section 12. 1. 11 we describe our first diffusion term in 
which cell-cell adhesions hinder motion, while in Section [2.1.21 we describe a different diffusion model where cell-cell 
adhesions promote motion. Section 12^ details our assumptions on MAPK activation in response to EGF treatment. 
In Section [2751 we then combine the models from Sections 12.1.ll and l2.1.21 with the observations on MAPK activation 
in Section [272[ to describe our models for EGF data in Sections |2.3.1l and l2.3.2L respectively. We then detail the entire 
model, including initial and boundary conditions, in Section [2.4! and outline our parameter estimation procedure in 
Section 12.51 

2.1. Models for diffusion 

In this section, we will derive two different models (Models H and P) for diffusion based off cell-cell adhesion. In 
order to do so, we first discretize our solution domain. We use a uniform grid for both the time and space domain, 
so the time interval may be denoted as t{j) = jAt, j = 0, ...,N — 1, where N denotes the total number of time 
points used, and x{i) = iAx,i = 0, ...,M — 1 where M denotes the number of spatial points used. We simplify 
notation by writing t{j) = tj, x(i) = Xi, and write our discretized solution as u(tj,Xi) = u^. 

2.1.1. Model H: cell-cell adhesions hinder migration 

For our first model on the effects of cell-cell adhesion on diffusion, we assume that cell-cell adhesions hinder 
migration through the drag of follower cells and accordingly denote it as Model H. In this scenario, the downregu- 
lation of cell-cell adhesions will promote migration. We create this model by modifying a previous cell migration 
study that incorporates cell-cell adhesion. 

In [l|, the authors developed a mathematical model that assumes cell-cell adhesions hinder cell movement. To 
derive the model, a probability transition is given for a cell density located at position Xi attempting to move 
forward to location by 


+ _ (1 - Ui+i){l - aui-i) 


( 2 . 2 ) 


where the first factor represents space filling (e.g., the cell density is more likely to move forward given a smaller 
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cell density in front of it). The parameter a denotes the (dimensionless) rate of cell-cell adhesion, and thus the 
second term in Equation represents our first hypothesis in which follower cells (denoted with Ui-i) hinder 

the movement of cells in front of them towards the wound. The transition probabilities , and are all 

defined analogously (with the minus signs denoting cells moving backwards). We note that Equation (|2.2p may be 
modified to include the influence of other aspects of cell migration, including cell-surface adhesion and shape change 

Q- 

Given these transition probabilities, Anguige and Schmeiser derive the continuum limit of the model by analyzing 
the change of cell density Ui over time as 


Ut = 7-+iUi_i -b T^+iU*+i - (t+ -b )w* (2.3) 

and then taking the limit as Ax —>■ 0. The derivation of the continuum limit of Equation (12.31) is presented in [ll 
and may be written as the dimensionless model 

Ut = ((1 -b 3a(u - 2 / 3 )^ - */3a)uo,)x, 


and we accordingly write the dimensional model as 


Ut = {{D -b 37 (u - 2 / 3 )^ - 'i/ 3 ^)ux)x- (Model H) 


(2.4) 


In Equation (12.4p . D and 7 denote the rates of diffusion and cell-cell adhesion with units microns^/^,. Hence the 
diffusion rate for Model H is 

fu =D + 37 (u - 2/3)2 - 4/3^. 


We note that backwards diffusion may occur with Equation (12.41) if 47 > 3D and address this inequality more after 
final model development in Section r2. 3. II 


2.1.2. Model P: cell-cell adhesions promote migration 

For our second model, we assume that cell-cell adhesions promote migration as leader cells pull the cells behind 
them forward and denote it as Model P. In this scenario, the upregulation of cell-cell adhesion will promote migration. 
We now denote the forward transition probability of cell density Ui as 


+ _ (1 - -b aui+i) 


(2.5) 
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where the first term again represent space filling, but the second term in Equation i2.5\) represents our second 
assumption in which leader cells (denoted with rti+i) promote movement towards the wound by pulling on the 
cells behind them. We note again that and are all defined analogously. Substituting the transition 

probability from Equation (12.51) into Equation (12.3p yields 


Ut 


(1 - Ui)(l + aui)ui-i + (1 - Ui){l + aui)ui+i - ((1 - Uz+i)(l + auj+i) + (1 - Uz-i)(l + auj- 

Aa;^ 

Ui_i + aul_^Ui - aui-iuf + Ui+i - 2ui - aufui+i + auiuf^^ 

Ax^ 

Ui —1 2u2 U 2 -I-I ^ rij ^Ui Ui — 


l))ui 


( 2 . 6 ) 


where we can recognize the first term on the right hand side as the standard central difference approximation to the 
second derivative. For the second term on the right hand side, we can simplify using Taylor series approximations 
as: Ri Ui + Ax ■ u' + ^^^(2 ■ u", iti_i th Ui — Ax ■ u' + ^^^2 ■ u", where primes denote a spatial derivative term. 

Substituting these terms into the second term on the right hand side of Equation (12.611 reveals 


Ui-i — 2ui + Ui+i I 2Ax‘^Ui{u'.j)‘^ + ufu'lAx^ + 0{Ax'^) 
A^2 +« 

= u" + a{2ui{ui)'^ + u^u'O + 0{Ax'^) 

= u'! + + 0{Ax^) 


resulting in the dimensionless continuum limit to be 


Ut = ((1 + au‘^)ua;)x, 


which we dimensionalize as 

Ut = {{D + ju‘^)ux)x, (Model?) (2.7) 

where D and 7 again denote the dimensionalized rates of cell diffusion and cell-cell adhesion, respectively. Thus 
the rate of diffusion for Model P is given as 

f:^=D+ju^. 


We note that a more general form of this equation has been used previously to model temperature fronts and can 
have traveling wave solutions under certain scenarios |4l| . 


2.2. MAPK activation and cell growth in response to EGF treatment 

While we have difficulty matching experimental data over the entire experiment, our models performed well in 
the time range of 10 hours after EGF treatment until the end of the experiment at t = 30 hours. In Figure I^TTl we 
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FRET Ratio and best fit iine 



Figure 2.1: Plotting the relative average FRET ratio of the cell population during the EGF experiment from t = 10 
to 30 hours, along with its best fit line, which is given by b{t) in the text. Note that this data has an value of 
0.98. 

have plotted the normalized average level of FRET ratio (indicative of MAPK activation in the population) during 
the experimental protocol. The activation levels after t = 10 appear linear, so we focus on matching data in this 
time interval (t = 10 — 30 hours) and use a linear term for MAPK activation. We have plotted a best-fit line to 
the normalized FRET ratio data from f = 10 — 30 hours in Figure im which is given by b{t) = —0.0675 -I- 0.0380t. 
Note that this data has an value of 0.98 to help justify this linear term. We thus set as a constant, c, in 
Equation (EH). Also note that b{t) > 1 after t « 28 hours, even though the data has been normalized, so we set 
m(t) = 1 for these values. 


2.3. Models for EGF data 

We now present the final models for collective migration in response to EGF treatment. We do so by assuming 
that the rate of cell-cell adhesion changes over time in response to the level of m activation. This requires a variable 
term, denoted by r(TO), instead of a constant value, 7 , for the rate of cell-cell adhesion. 

2.3.1. Model H: MAPK stimulates collective migration through decreases in drag strength 

To conclude Model H, we now extend fff in Equation (12.41) to include changes in cell-cell adhesion in response 
to different levels of MAPK activation. This model is consistent with the hypothesis that MAPK activity stimulates 
collective migration through decreases in the drag strength of follower cells. We use a linear cell-cell adhesion term 
as a function of MAPK activity, m, by F(m) = 7(1 — m), for some cell-cell adhesion constant 7 > 0. Note that 
0 < m < 1 so that r(TO) > 0 . 

Accordingly, our final description for Model H in response to EGF treatment is given by 


Ut{t, x) = {{D 3F(m)(u — 2 / 3)2 _ '^/3V{m))ux)x^ ixit^t) = c (Model H for EGF) 


(2.8) 





where r(TO) is given above. We note that if AD < 3r(m), the above equation will exhibit backwards diffusion, 
which is not well-defined. Fortunately, all parameter sets in this study do not satisfy this inequality for the entirety 
of the simulation. 


2.3.2. Model P: MAPK activation stimulates collective migration by increases in pulling strength 

We extend from Equation (|2.7I) to include changes in cell-cell adhesion in response to different levels of MAPK 
activation for Model P in response to EGF data. In this scenario, we assume that MAPK activity stimulates collective 
migration through increases in the pulling strength of leader cells. Accordingly, we model the cell-cell adhesion level 
as a function of m by r(TO) = 7 m, for some constant 7 . As a result. Model P for EGF-treated cells is given by 

Ut{t,x) = ({D+ T{m)u^)ux)x + ku{^ — u), mt(t) = c, (Model P for EGF) (2.9) 

where P)™) is given above. 


2 . 4 . Complete models, initial and boundary conditions 

We’ve now introduced the two models under consideration in this study. We briefly review the two models, 
along with their terms and assumptions in Table [T] A list of all parameters used throughout this study is also given 
in Tabled] 


Model features 

H 

P 

Cell-cell adhesions 

Hinder migration 

Promote migration 

Simulation 

mock 

EGF 

mock 

EGF 

Equation 

(12.41) 

(12.81) 

(12.71) 

(12.91) 

MAPK impact on 
cell-cell adhesions 

n/a 

Downregulation 

n/a 

Upregulation 

r(m) 

7 

7(1 — m) 

7 

7 m 


Table 1: Summary of equations and assumptions used for Models H and P. 


If we let w denote the location of the wound, then we would ideally use an initial condition that is nonzero for 
all a: < re to imitate the presence of cells and zero for all x > w to imitate the wound. Scratching the cell population 
with the pipette tip, however, initially shocks the cells and causes slow initial movement. For this reason, we use a 
nonzero initial condition for all x < u; = re — 180 microns and a zero initial condition for all x > w. This alteration 
to the initial condition allows us to match the wave phenomena of both the model and data. We compute the 
nonzero portion of the initial condition by interpolating the data profile at the initial time point. Thus, our initial 
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Parameter 

Description 

(units) 

Model H 

Model P 

Determined 

by 

Mock 

EGF 

Mock 

EGF 

D 

Baseline rate of diffusion 

(microns’^/hr) 

19,770 

177,940 

2.19 

0.14 

Fitting 

c 

Rate of MAPK activation 
(hr-i) 

0 

0.038 

0 

.038 

h{t) 

7 

Rate of adhesion between adjacent cells 
(microns^/hr) for mock 
(microns'/hr.m) for EGF 

1.68 

970 

30,614 

535,810 

Fitting 

w 

Location of wound 
(microns) 

1550 

1620 

1550 

1620 

From data 


Table 2: List of parameters for the models, along with the parameter values used and how we obtained these values. 


X < w 
otherwise 

( 2 . 10 ) 

where to is our starting time for simulations (e.g., to = 10 hours) and ip denotes the interpolated, normalized initial 
data profile values at time to. In Figure we’ve included some typical video snapshots at t = 7.5, 20 hours, along 
with their resulting density prohle. In order to obtain this population prohle from data, we sum over the vertical 
axis for each image matrix and then normalize over each time step. The linear interpolation of this initial profile (at 
t = to) is our '(/'(x) function used for the initial condition. Because of the linear MAPK activation phase between 
t = 10 — 30 hours, we will simulate the model from t = 7.5 to 30 hours, where we also run the model for an extra 
2.5 hours due to its initial high speeds. We also include the leading edge locations of the two experiments in the 
bottom portion of Figure [221 whose calculation is described in the following section. 

We use zero Neumann boundary conditions to simulate no flux conditions at the walls of the well plate. Note 
that in experimental videos, we only observe a 3.24 mm x 3.24 mm field of view, whereas the experimental domain is 
actually 7 mm long and 5 mm wide, so there is plenty of space for the cells to move in either direction. Accordingly, 
we extend our computational domain 1.62 mm behind the field of view and 3.24 mm in front of the field of view to 
ensure the boundaries do not affect the model simulations. 


condition is given by 


u{to,x) = 


1 p{x)] 


m{to) = b{to). 


2.5. Parameter estimation 

In order to fit our model simulations to data on wound healing, we will use an inverse problem procedure to 
estimate the parameter vector q = [£>, a]^ for both models for both mock and EGF data. We do so by comparing 
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Location fmicrons) 


(a) Experimental snap shots of the experiment and resulting 
profiles for the mock experiments. Top row: Images of the 
wound healing experiment at f = 7.5 and 20 hours. Bottom 
row: Resulting profiles for the cell population in the image. 


t = 7.5 hours 


t = 20 hours 





g 0.6 

3 0.4 

E 
w 


1500 3000 

Location imicronsl 




1500 3000 

Location imicronsl 


(b) Experimental snap shots of the experiment and resulting 
profiles for the EGF experiments. Top row: Images of the 
wound healing experiment at t = 7.5 and 20 hours. Bottom 
row: Resulting profiles for the cell population in the image. 


Leading edge locations for mock and EGF data 



(c) Locations of the mock (blue) and EGF (red) experimental 
leading edges over time using Equation 112. lil t , which identifies 
the first point where u(tj,Xi) < 0.3 for each time point tj, as 
is described in Section 12.51 Note that the mock leading edge 
data may take the same value for several consecutive time 
points because it moves slowly, and the camera discretizes 
the experimental domain. 


Figure 2.2 


11 





































the locations of the leadin g edge f or both the experimental data and model simulations, as has been done in previous 


studies on wound healing 


25 


30|. In this work, we will define the leading edge as the location in which the relative 
density is equal to a certain value, although we note that it has different definitions elsewhere As seen in Figures 
3.31 and r3.4l a density of 0.3 can locate the leading edge well. Accordingly, we define imfi.sit, <?) as the leading edge 
location over time for the model using q as an estimate for the true q, denoted by 


^m,o. 3 (i, q) = {x \u{x, t) = 0.3 , q=q} . 


Note that with numerical computations, we cannot precisely estimate the leading edge in this manner. Accordingly, 
we estimate the leading edge to be the first location where the relative density is below the set threshold. For 
instance, if q 3 (g) denotes the leading edge location at time then we approximate it as: 


C.o.3('7) ~ C,0.3(9) = {xi K < 0.3 and > 0.3 for all /c < i , q=q} . ( 2 . 11 ) 

We use analogous definition for the estimate of the leading edge of the normalized data over time, which is denoted 

^z>,o.3- 

As a means to estimate q using our models, we will implement an inverse problem in which we minimize a cost 
function . In this work, we will use the cost function given by an ordinary least squares estimate: 

N 

^(9l = ElC.o.3(9)-^\o.3p, 

n—1 

where N is the number of time points considered. To find the value of q that minimizes the cost function, we use 
the Nelder-Mead algorithm as implemented in MATLAB’s fminsearch command. 


3. Results 

To investigate the performance of Models H and P in describing the wound healing process, we now compare 
both models to mock and EGF experimental data. To compare these prediction to the experimentally determined 
results, we use images of the cell populations over time, such as those depicted in Figure 12.21 We use several 
different criteria for comparison to the experimental data. In Section 13.11 we compare the leading edge locations of 
the model simulations and experimental data for both of the mock and EGF cases. In Section [321 we compare the 
wave speed values of the EGF model simulations to the experimental data wave speed values. In Section 13.31 we 
directly compare the experimental snapshots to the model profiles over time. We conclude the results section by 
using Model P to investigate the time to wound closure as a function of the wound area in Section 13.41 
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3.1. Leading edge fitting 

Using the protocol described in Section 12.51 we fit both models to both mock and EGF data. The resulting 
leading edge simulations are depicted in Figure 13.11 and the parameters used to obtain these plots are given in 
Table [2] The red dots in Figure 13.11 represent the location of the experimental leading edge over time, while the 
green and blue lines represent the leading edge locations for Models H and P, respectively. The thin black dash-dot 
lines represent one standard deviation of the data, and we describe their computation in the appendix. Note that 
the scale for the two images is very different, as the leading edge in FGF simulations moves much faster and farther 
than the leading edge in mock simulations. Recall that mock leading edge data sometimes takes the same value for 
several consecutive time points because it moves slowly, and the camera discretizes the experimental domain. 

In the mock simulations (depicted in the left image of Fieure lXTl) . both models appear to match the leading edge 
data in a similar manner. In particular, after t = 18 hours, the two leading edge simulations appear superimposed, 
suggesting that these two very different models can yield similar results for predicting the the leading edge location. 
Table [5] summarizes, however, that these two models use very different simulations to obtain these similar results. 
Note that Model H relies primarily on diffusion for movement, as it uses a large baseline diffusion constant of 
D = 19, 770 microns^/hr, compared to Model P, which uses a much smaller value oi D = 2.19 microns^/hr. On 
the other hand. Model P uses a large rate of cell-cell adhesion, as it estimates 7 = 30, 614 microns^/hr while Model 
H estimates 7 = 1.68 microns^/hr for its rate of cell-cell adhesion. Thus we see that Model H uses a high rate of 
diffusion to match the mock data, whereas Model P uses a high rate of cell-cell adhesion to match this data. These 
results are not surprising based on the derivations of these two models, as cell-cell adhesions hinder movement for 
Model H but promote movement for Model P. 

Extending the same parameter estimation methodology to EGF data, we see on the right hand side of Figure 
o that Model P fits the experimental data better than Model H. Model P matches the experimental data well at 
the initial time point of t = 10 hours and continues to do so until the hnal time point. On the other hand, Model H 
initially overestimates the leading edge from t = 10 — 20 hours (so much so that it even fails to match the data to 
within one standard deviation at times) and then underestimates the leading edge location from t = 20 — 30 hours. 
We thus observe that Model H is unable to accurately simulate the EGF data, suggesting that MAPK activation 
may not decrease drag strength during wound healing when the cell population is treated with EGF. 

Looking at the parameter values for Model P simulations in Table [5J we see that the estimate for 7 increases 
significantly in response to EGF treatment. More surprising, however, is the decrease in the rate of baseline cell 
diffusion that we observe between the mock and EGF simulations for Model P, from D = 2.19 microns^/hr in the 
mock experiment to I? = 0.14 microns^/hr in the EGF experiment. Due to the relative insensitivity of Model P 
to the parameter D (results not shown), we can likely consider these values of D as equivalent. Hence, MAPK 
activation appears to have a negligible effect on the baseline rate of diffusion. This is an interesting observation, as 
one may initially suspect that increases in MAPK activation increase the baseline rate of diffusion as a means to 
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Mock Leading Edge Simulations 


EGF Leading Edge Simulations 




Figure 3.1: Leading edge simulations for both mock- (left) and EGF- (right) treated cell populations. The red dots 
represent the experimental leading edge locations over time, and the solid green and dashed blue lines represent 
the leading edge locations for Models H and P, respectively. Note that while both models can match the mock 
data relatively well, model P appears to fit EGF data better than model H. The black dash-dot lines denote one 
standard deviation of the leading edge locations, and their computation is described in the appendix. 

stimulate migration, as we observe for Model H, whose estimates of D change from 19,770 to 177,940 microns/hr^ 
in response to EGF treatment. 


3.2. Wave speed comparisons 

A previous study on predicting the time for wound closure in intestinal epithelial cells used instantaneous velocity 
of the wound edge to compare model simulations with experimental data In a similar fashion, we compare wave 
speed simulations of both the experimental data and the best-fit model simulations from the previous section. We 
use a forward difference calculation on the leading edge locations over time to estimate the model simulation wave 
speeds, w"!^ for various normalized cellular density levels (given by /3 = 0.2, 0.3, 0.4): 


w 


n _ 

m,(3 


pn+l _pn 
^m,f3 ^rn,0 


At 


and the data wave speeds are calculated analogously. In this equation, ^ denotes the 4-point moving average 
(computed using MATLAB’s smooth function) of the leading edge location at time tn using /3 as the leading edge 
threshold. We use ^ instead of ^ in this calculation to obtain realistic wave speed values. Note that the data 
has a very sharp front, which causes similar wave speed values for all densities considered from t = 10 — 20 hours, 
which is demonstrated in Figure l3.2bl 

The results for EGF scenario model and data wave speeds for /3 = 0.2, 0.3,0.4 are depicted in Figure I5.2al in 
green, blue, and red, respectively. Note that both model wave speeds appear to have two different behavior phases: 
an initial fast phase (likely due to the presence of empty space), followed by a sustained, slowly decreasing phase. 
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The previous study also observed these two phases of behavior and sometimes also observed an increase in speed 
towards the end of the simulation due to wound closure. Compare these changes in phase behavior with that of 
the the data, which initially has a constant velocity but begins to slow down after t = 20 hours. The initial high 
velocity of the model simulations (and lack thereof in the experimental data) is why we let the model simulations 
run for 2.5 hours before comparing their leading edges with the experiments. 

In comparing the two models’ wave speed simulations to the data, we see that Model H consistently underes¬ 
timates the wave speed for /3 = 0.3, 0.4, and overestimates the wave speed for /3 = 0.2. On the other hand. Model 
P can match the faster phase of the data for /3 = 0.2 (around 55 microns/hr) from t = 17 — 30 hours, and the 
slower phases of the data (between 35 and 45 microns/hr) for /3 = 0.3 and /3 = 0.4 from f = 20 — 29 hours. Overall, 
Model P tends to match the data better than Model H for all time points. In Figure I3.2b[ we superimpose all 
data and model wave speed values to obtain a better picture of the overall population wave speeds. Here we again 
observe more agreement between the experimental data and Model P than we observe for Model H. From these 
simulations, it is not surprising that Model P can more accurately fit the leading edge data than Model H, as its 
best-fit simulation better matches the wave speeds of the data. 

3.3. Comparison to data snap shots 

We next compare the model profiles to experimental snapshots over time for both mock and EGF simulations. 
These images are depicted in Figures 13751 and lX^ for mock and EGF data, respectively, for t = 10,17.5, 25, and 29.5 
hours. In each snapshot, the top frame represents the experimental image of the cell sheet during the wound healing 
experiment, and the bottom frame demonstrates the model simulations for both Models H and P. In each snapshot, 
the experimental leading edge is depicted with a black line in the top frame and a black dot in the bottom frame. 

In Figure [3731 we note that both Models H and P match the experimental leading edge well, as we would expect 
from Figure [37T] However, the two resulting profiles are very different, as Model H has a gradually decreasing profile 
in comparison to the steep population front of Model P. Gonsidering our previous observation that Model H uses 
a high rate of diffusion to match the leading edge, whereas Model H uses a high rate of cell-cell adhesion, these 
two different model profiles are not surprising. Note that this gradually descending cell profile causes Model H to 
predict a significant presence of cells in locations that are clearly empty in the experiment. For example, at t = 29.5 
hours. Model H predicts a cellular density of about u{t = 29.5, a; = 2500) = 0.25 (one quarter of the maximum 
cellular density at a; = 2500 microns), whereas the experimental cell sheet has clearly not yet reached this location, 
and its leading edge is approximately 500 microns further back. We believe this difference in profile front descent 
allows Model P to better predict the leading edge location and various wave speeds of the experimental wound 
healing data. 

Figure [37T] depicts the same figure for the EGF experiments and simulations. We again notice that both profiles 
tend to agree near the leading edge location, but yield very different profile simulations. Model H’s profile again 
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Leading edge wavespeeds, beta = 0.2 



10 15 20 25 

Time (hours) 


Leading edge wavespeeds, beta = 0.3 



Leading edge wavespeeds, beta = 0.4 



(a) Wave speed comparisons for EGF scenarios for Model H (solid line), Model P (dashed line), and experimental data (various shapes). 
The wavespeeds for (3 = 0.2 are compared in green, for (3 = 0.3 are in blue, and for /3 = 0.4 are in red. Note that model P tends to 
match the data better than Model H, and that Model H consistently over and underestimates almost all wave speed values from the 
data for /3 = 0.2 and /3 = 0.4, respectively. 



(b) Superposition of all density wave speeds with data. Model H is depicted on the left hand side, and Model P is depicted 
on the right hand side. Note that Model P agrees more with the data than Model H. 


Figure 3.2: Wave speed simulations for the two models in comparison with data. 


16 



















Mock Model comparison, Time = 10 hours 



Mock Model comparison, Time = 17.5 hours 




Mock Model comparison. Time = 25 hours 



Mock Model comparison, Time = 29.5 hours 





Figure 3.3: Comparisons of experimental images (upper frames) to model profiles (lower frames) for the mock 
experiments over various time snapshots. The black line in the upper frames and dot in the bottom frames depict 
the experimental leading edge over time. Note the gradually decreasing cell front of Model H in comparison to the 
sharp front of Model P. 


17 














































EGF Model comparison, Time = 10 hours 
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EGF Model comparison, Time = 17.5 hours 
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EGF Model comparison. Time = 25 hours 
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Figure 3.4: Comparisons of experimental images (upper frames) to model profiles (lower frames) for the EGF 
experiments over various time snapshots. The black line in the upper frames and dot in the bottom frames depict 
the experimental leading edge over time. Note the gradually decreasing cell front of Model H in comparison to the 
sharp front of Model P. Observe how this gradual population frony causes Model H to predict a nonzero value for 
u{t = 17.5, X = 3000), even though there are no cells present at this time and location in the experiment. 


descends much more gradually than the steep front of Model P. Similarly, Model H also predicts the presence of 
cells in areas that are clearly empty in the experimental snapshot, such as at a; = 3000 microns at t = 17.5 hours. 
We note, however, that Model P also predicts the presence of cells in clearly empty areas from the EGF experiment 
(such as past x = 3000 microns at t = 25 hours). This suggests that more sophisticated modeling endeavors are 
still needed in the future to better capture the leading edge dynamics and match the experimental profiles of wound 
healing data. Overall, however, Model P can accurately match the leading edge dynamics well for both mock and 
EGF simulations, whereas Model H is considerably flawed in several aspects. These results support our second 
hypothesis that MAPK activation in response to EGF treatment stimulates collective migration by increasing the 
pulling strength of leader cells. 
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Figure 3.5: Using Model P to predict the time to wound closure as a function of wound area. The blue diamonds 
denote closure time for mock cells, whereas the red dots denote closure time for EGF-treated cells. Note that the 
time to wound closure for mock-treated cells appears quadratic, which matches results from . The time to wound 
closure appears to become linear in response to EGF treatment. 
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3.4- Time to wound closure 

A common theme in wound healing studies is to investigate how different wound aspects, including size, shape, 
and aspect ratio, affect the time to wound closure 2, Il5| . Times for wound closure as a function of wound area 




were analyzed in [2] for different wound shapes and sizes. In a similar fashion, we can use Model P to predict the 
closure times for various different wound sizes, as is depicted in Figure 13.51 Because this model is simulated in 
one dimension, we only use it on rectangular wound shapes, where we calculate the wound area by multiplying the 
length of the wound by the width of the well (4 microns). We also define the time to wound closure as the time 
it takes for the leading edge to reach the end of the field of view (3,240 microns), or such that 3 = 240 

microns. Notice in this figure how the time to wound closure for mock cells appears approximately quadratic, which 
would be consistent with j^, in which the authors suggested that the rate of change of wound area is proportional 
to the area of the wound in lEC cells. In response to EGF treatment, however, we see this curve change to a much 
faster and apparently linear term. We can validate these predicted times in future studies by experimenting with 
wounds of various sizes. We may also investigate nonrectangular wounds in the future by extending our model to 
two dimensions. 


4. Conclusions and Discussion 

In this study, we have developed two mathematical models representing the spread of a cellular sheet during 
wound healing in response to EGF treatment. These two models are used to test two hypotheses on the effects of 
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MAPK activity on cell-cell adhesion during collective migration. Model H assumes that MAPK activity stimulates 
collective migration through decreases in the drag strength of follower cells, whereas Model P assumes that MAPK 
activity stimulates collective migration through increases in the pulling strength of leader cells. Model P matches 
several aspects of the experimental data better than Model H, ultimately suggesting the validity of our second 
hypothesis. From our resulting parameter values, we also observe that the main effect of MAPK activity on 
cell migration is an increase in the rate of cell-cell adhesion, and that it has a negligible effect on the baseline 
rate of diffusion. As one may initially suspect that increases in MAPK activation stimulate the baseline rate of 
diffusion to promote the rate of migration, we plan to further investigate this observation in future experimental 
and computational studies due to its nonintuitive nature 

These conclusions have implications in our understanding of the EMT, in which cells detach from an epithelial cell 


population as mesenchymal cells 


flEl. 


Cells undergoing EMT tend to gain migratory and invasive properties, and 


n 


hence the EMT is a driving factor in both normal embryological development and cancer metastasis [7 


221. There 


are two main types of EMT: incomplete and complete. During incomplete EMT, cell-cell junctions remain intact 
between certain invasive cells, and collective strand-like migration is observed. Complete EMT is characterized 
by the loss of most cell connections, and individual cells detach from the population to begin their own migration 
pattern. Hence, in either form of EMT, we see that cells detach from the population and can be considered smaller 
cell subpopulations (Composed of one or several cells), effectively increasing the prevalence of leader cells. Since cell 
pulling appears to be the dominant mechanism of collective migration, we suggest that HaCaT cells in response to 
EGE treatment may express incomplete EMT behaviors so that an increased number of leader cells can effectively 
guide smaller cell subpopulations via pulling. We may use this observation to better understand incomplete EMT 
dynamics through the experimental system and mathematical analysis, along with an understanding of any possible 
connections between incomplete and complete EMT. For example, if any sort of chemical treatment or oncogenic 
perturbations will convert the system from incomplete to complete EMT or vice versa. 


In Figure 14.11 we have depicted several data and model contours over time. From these contours, we see that 
Model P is able to accurately match the leading edge data for both mock and EGF data, but it is still far from 
matching the entire data profile accurately over time. Future endeavors will aim to match the entire population 
profile instead of just the leading edge. We note that the role that the MAPK signaling cascade on cell-surface 
adhesion during collective migration is another important question that requires investigation. Further exploration 
into this may aid in better matching the entire cell population profile. We also note that we used a simple linear 
term for average MAPK activation levels in our model. The slope of this line very likely depends on the amount of 
EGF treatment in the population, however, so future studies will include more sophisticated modeling terms that 
may involve EGF ligand, its intracellular precursor, and reactive oxygen species, as has been done in a previous 


model on MAPK activation during wound healing 


38j . EGF ligand secretion by leader cells likely aids in collective 
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Figure 4.1: Contours for the two models (top two rows) against data (bottom row). The x—axis depicts the spatial 
location and the ?/—axis depicts time to demonstrate how the profiles change over time for both model simulations 
and the actual experiment. The colorbar depicts the cellular density levels (i.e., blue denotes a low density and red 
denotes a high density). 


migration as ligand diffusion to follower cells will further promoting cell pulling in the population. 

In this work, we have developed a simple model to investigate how EGF treatment influences collective migration 
during wound healing in HaCaT cells. We demonstrate that cell-cell adhesion plays a critical and activating role in 
collective migration through leader cell pulling. This model may have implications in understanding EMT dynamics, 
including how cells may transition from incomplete to complete EMT and vice versa. 
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Appendix 

Numerical implementation 

Note that for numerical implementation, we use a method of lines approach with MATLAB’s odel5s function 
to simulate the entire model. For spatial discretization, we use the second order scheme for convection-diffusion 
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Figure 4.2: Splitting up the y-axis into twenty bins and obtaining several estimates for the leading edge. On the 
left hand image, we’ve demonstrated how we split up the axis to calculate several different leading edge locations 
(denoted by black vertical lines). On the right hand image, we show the histogram of the leading edge locations of 
the left hand image . 


equations (without convection in our case) from 


28| given by 


• u\ _ Pj-y2{t) 


where (t) is an approximation to the diffusive flux, given by 


where Q{u,Ux) 
simulations. 


Pj+^/2{t) — 2 


Q ( Uj{i)^ 


Uj+i{t) - Uj{t) 

Ax 


denotes the cellular diffusion rate. Simulations 


not shown demonstrate quick convergence in our 


Computing standard deviations 

In order to compute the standard deviation of the location of the leading edge at each time point, we split 
the 1 /—axis of each image into R bins, and compute the leading edge for each of the separate bins to obtain R 
estimates of the leading edge location in the population. In Figure W/H for example, we’ve split the y—axis into 
twenty separate bins on the right hand side and calculated the leading edge of the data for each one. Hence, for 
each time point during the experiment, we compute the sample mean of the leading edge locations (call it /t) then 
compute the sample standard deviation to obtain one standard deviation. On the right side of Figure 14.21 we’ve 
also depicted an example histogram for the leading edge locations when using twenty bins. 
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